function [alphas,betas,R2s, mats_yrs] = FBregs_real_ver2(gesol, Phi_Nz, stationary_prob_Nz)
% regression: (n-1)*(y(n-1,t+12)-y(n-1,t))=a+b*fwd_spread(n,t)+err

    mats_yrs = 2:5;
    
    alphas = zeros(size(mats_yrs));
    betas = zeros(size(mats_yrs));
    R2s = zeros(size(mats_yrs));
    
    for i = 1:numel(mats_yrs)
        
        mat_yr = mats_yrs(i);
        mat_mth = mat_yr*12;
        
        % RHS: forward spread = p(n-1) - p(n) + p(1)
        idx_1y = find(gesol.real_bonds.maturities==12);
        idx_ny = find(gesol.real_bonds.maturities==12*mat_yr);
        idx_nless1y = find(gesol.real_bonds.maturities==12*(mat_yr-1));
        
        if isempty(idx_1y) || isempty(idx_ny) || isempty(idx_nless1y)
            error('Some yields are not available');
        end
        
        rhs_Nz = log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y))) ...
            - log(squeeze(gesol.real_bonds.prices(:,:,:,idx_ny))) ...
            + log(squeeze(gesol.real_bonds.prices(:,:,:,idx_1y)));
        
        % LHS = (n-1)*y(n-1,t+1) - (n-1)*y(n-1,t)
        % LHS1 = -(n-1)*y(n-1,t), LHS2 at h=12 is (n-1)*y(n-1,t+1)
        
        lhs1_Nz = log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y)));
        lhs2_Nz = -log(squeeze(gesol.real_bonds.prices(:,:,:,idx_nless1y)));
        
        mean_rhs = sum(stationary_prob_Nz(:).*rhs_Nz(:));
        mean_lhs1 = sum(stationary_prob_Nz(:).*lhs1_Nz(:));
        mean_lhs2 = sum(stationary_prob_Nz(:).*lhs2_Nz(:));
        var_rhs = sum(stationary_prob_Nz(:).*(rhs_Nz(:).^2)) - mean_rhs^2;
        var_lhs1 = sum(stationary_prob_Nz(:).*(lhs1_Nz(:).^2)) - mean_lhs1^2;
        var_lhs2 = sum(stationary_prob_Nz(:).*(lhs2_Nz(:).^2)) - mean_lhs2^2;
        cov_lhs1_lhs2 = sum(stationary_prob_Nz(:).*lhs1_Nz(:).*((Phi_Nz^12)*lhs2_Nz(:))) - mean_lhs1*mean_lhs2;
        var_lhs = var_lhs1 + var_lhs2 + 2*cov_lhs1_lhs2;
        cov_rhs_lhs1 = sum(stationary_prob_Nz(:).*rhs_Nz(:).*lhs1_Nz(:)) - mean_rhs*mean_lhs1;
        cov_rhs_lhs2 = sum(stationary_prob_Nz(:).*rhs_Nz(:).*((Phi_Nz^12)*lhs2_Nz(:))) - mean_rhs*mean_lhs2;
        
        betas(i) = (cov_rhs_lhs1 + cov_rhs_lhs2)/var_rhs;
        alphas(i) = mean_lhs1 + mean_lhs2 - betas(i)*mean_rhs;
        R2s(i) = (betas(i)^2)*var_rhs/var_lhs;
        
    end
    

end